/*
  Copyright 2006-2011 Patrik Jonsson, sunrise@familjenjonsson.org

  This file is part of Sunrise.

  Sunrise is free software; you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation; either version 3 of the License, or
  (at your option) any later version.

  Sunrise is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with Sunrise.  If not, see <http://www.gnu.org/licenses/>.
  
*/

// $Id$

#include<iostream>
#include "CCfits/CCfits"
#include "counter.h"
#include "polychromatic_grid.h"
#include "optical.h"
#include "scatterer.h"
#include "constants.h"
#include "dust_model.h"
#include "emission.h"
#include "emission-fits.h"
#include "full_sed_emergence.h"
#include "full_sed_grid.h"
#include "emergence-fits.h"
#include "mcrx-units.h"
#include "grid_factory.h"
#include "grid-fits.h"
#include "shoot.h"
#include "fits-utilities.h"
#include "boost/lexical_cast.hpp"
#include "boost/shared_ptr.hpp"
#include "density_generator.h"
#include "preferences.h"
#include "model.h"
#include "grain_model_hack.h"
#include "p04_grain_model.h"
#include "wd01_grain_model.h"
#include "xfer.h"
#include "imops.h"
#include "ir_grid.h"
#include "equilibrium.h"
#include "boost/lambda/lambda.hpp"
#include "boost/program_options.hpp"

//using namespace blitz;
using namespace mcrx;
using namespace CCfits;
using namespace std; 
using namespace boost::lambda;
using boost::lexical_cast;
namespace po = boost::program_options;

typedef local_random T_rng_policy;
typedef scatterer<polychromatic_scatterer_policy, T_rng_policy> T_scatterer;
typedef T_scatterer::T_lambda T_lambda;
typedef T_scatterer::T_biaser T_biaser;
typedef dust_model<T_scatterer, cumulative_sampling, T_rng_policy> T_dust_model;
typedef full_sed_emergence T_emergence;

typedef full_sed_grid<adaptive_grid> T_grid;
typedef ir_grid<adaptive_grid, cumulative_sampling, local_random> T_ir_grid;


class pascucci_factory: public grid_factory<T_grid::T_data> {
public:
  const T_densities rho0_, n_tol_;
  const T_float ri, ro, rd, zd, r_tol_;
  const int maxl, n_lambda;
  
  pascucci_factory (const T_densities& rho0, const T_densities& n_tol,
		    T_float r_tol,
		    int ml, int nl):
    ri(1), ro(1000), rd (ro/2), zd(ro/8),
    rho0_(rho0), n_tol_(n_tol), r_tol_(r_tol),
    maxl(ml), n_lambda (nl) {};

  T_densities calc_rho(const vec3d& p) {
    T_densities rho;
    resize_like(rho, rho0_);
    rho = 0;
    const T_float r = sqrt(p[0]*p[0] + p[1]*p[1]); 
    const T_float z = p[2];
    if ((r < ro) && (r > ri)) {
      // we are in the allowed range of radii, calculate the density
      const T_float h=zd*pow(r/rd,1.125);
      const T_float f2=exp(-constants::pi/4*pow(z/h,2));
      rho = rho0_ *(rd/r)*f2;
    }
    return rho;
  };

  virtual bool refine_cell_p (const T_cell& c, int level) {
    if (level>maxl) 
      return false;
    //if (level<=4) return true; // a minimum amount of "prerefinement"
    // refine until cell column density < n_tol

    // simone uses maxlevel=10+log(r/rmin)
    const vec3d min = c.getmin();
    const vec3d max = c.getmax();
    /*    
    T_densities rhoav
      ((calc_rho(c.getcenter())+
	calc_rho(min)+
	calc_rho(vec3d(min[0],min[1],max[2]))+
	calc_rho(vec3d(min[0],max[1],min[2]))+
	calc_rho(vec3d(min[0],max[1],max[2]))+
	calc_rho(vec3d(max[0],min[1],min[2]))+
	calc_rho(vec3d(max[0],min[1],max[2]))+
	calc_rho(vec3d(max[0],max[1],min[2]))+
	calc_rho(max))/9);
    */
    const T_float sz=sqrt(dot(c.getsize(),c.getsize()));
    T_densities rhoav
      (-log((exp(-sz*calc_rho(c.getcenter()))+
	     exp(-sz*calc_rho(min))+
	     exp(-sz*calc_rho(vec3d(min[0],min[1],max[2])))+
	     exp(-sz*calc_rho(vec3d(min[0],max[1],min[2])))+
	     exp(-sz*calc_rho(vec3d(min[0],max[1],max[2])))+
	     exp(-sz*calc_rho(vec3d(max[0],min[1],min[2])))+
	     exp(-sz*calc_rho(vec3d(max[0],min[1],max[2])))+
	     exp(-sz*calc_rho(vec3d(max[0],max[1],min[2])))+
	     exp(-sz*calc_rho(max)))/9)/sz) ;
    
    // refine cell if cellsize/distance from source > r_tol_ or if it
    // is within one cellsize of the boundary (assuming some of the
    // densities are nonzero)
    vec3d cylr=c.getcenter();
    cylr[2]=0;
    if ((any(rhoav)>0) && 
	//( (dot(c.getsize(), c.getsize()/dot(c.getcenter(), c.getcenter())) >
	( (dot(c.getsize(), c.getsize()/dot(cylr,cylr)) >
	   r_tol_*r_tol_) ||
	  (2*abs(sqrt(dot(cylr,cylr)-ri))<
	   sqrt(dot(c.getsize(),c.getsize())))
	  ))
      return true;

    // refine cell if column density exceeds n_tol
    // length required for tau 1
    //T_densities abslen(64*n_tol_/rhosum);
    //if (condition_all(abslen>50))
      return condition_all (rhoav*rhoav*dot(c.getsize(),c.getsize()) > 
			    n_tol_*n_tol_);
      //else
      // hopelessly optically thick, just go by the radial criterion
      //return false;
  };
  virtual bool unify_cell_p (T_grid::const_local_iterator begin,
			     T_grid::const_local_iterator end) {
    return false;
    /*
    vec3d center = begin->getcenter();
    T_float vol = begin->volume();
    T_densities s ((begin++)->data()->get_absorber().opacity());
    T_densities ssq (s*s);
    int n=1;
    while (begin != end) {
      center += begin->getcenter();
      vol += begin->volume();
      T_densities o = (begin++)->data()->get_absorber().opacity();
      s+=o;
      ssq+=o*o;
      n++;
    }

    const T_densities mean (s/n);
    const T_densities stdev2 (ssq/n- mean*mean);
    const T_float r2=dot(center,center)/(n*n);
    return condition_all ((mean==0) || ((stdev2/(mean*mean))<tol) || 
			  (stdev2<r2*r2/(vol*vol)*tabs*tabs));
    */
  };

  virtual T_data get_data (const T_cell& c) {
    const vec3d min = c.getmin();
    const vec3d max = c.getmax();
    /*
    T_densities rhoav
      ((calc_rho(c.getcenter())+
	calc_rho(min)+
	calc_rho(vec3d(min[0],min[1],max[2]))+
	calc_rho(vec3d(min[0],max[1],min[2]))+
	calc_rho(vec3d(min[0],max[1],max[2]))+
	calc_rho(vec3d(max[0],min[1],min[2]))+
	calc_rho(vec3d(max[0],min[1],max[2]))+
	calc_rho(vec3d(max[0],max[1],min[2]))+
	calc_rho(max))/9);
    */
    const T_float sz=sqrt(dot(c.getsize(),c.getsize()));
    T_densities rhoav
      (-log((exp(-sz*calc_rho(c.getcenter()))+
	     exp(-sz*calc_rho(min))+
	     exp(-sz*calc_rho(vec3d(min[0],min[1],max[2])))+
	     exp(-sz*calc_rho(vec3d(min[0],max[1],min[2])))+
	     exp(-sz*calc_rho(vec3d(min[0],max[1],max[2])))+
	     exp(-sz*calc_rho(vec3d(max[0],min[1],min[2])))+
	     exp(-sz*calc_rho(vec3d(max[0],min[1],max[2])))+
	     exp(-sz*calc_rho(vec3d(max[0],max[1],min[2])))+
	     exp(-sz*calc_rho(max)))/9)/sz) ;
    //const T_densities rhoav = calc_rho(c.getcenter());
    return T_data ( rhoav, vec3d(0,0,0), T_lambda (n_lambda) );
  }
    
  virtual int n_threads () {return 1;};
  virtual int work_chunk_levels () {return 5;};
  virtual int estimated_levels_needed () {return 0;};
};


class uniform_factory: public grid_factory<T_grid::T_data> {
public:
  const T_densities rho_, n_tol_;
  const T_float rmax_;
  const int n_lambda;

  T_densities get_rho(const T_cell& c) {
    // return 0 if outside of radius rmax
    T_densities rho = (dot(c.getcenter(), c.getcenter())< rmax_*rmax_) ? rho_ : T_densities(rho_*0.);
    return rho;
  };

  uniform_factory (const T_densities& rho, T_densities n_tol, 
		   int nl, T_float rmax):
    rho_ (rho), n_lambda (nl), n_tol_(n_tol), rmax_(rmax) {};
  
  virtual bool refine_cell_p (const T_cell& c, int level) {
    // cell is refined to keep column density below n_tol
    T_densities rho = get_rho(c);
    T_densities n2 (dot(c.getsize(), c.getsize())*rho*rho);
    bool refine = all( n2 > (n_tol_*n_tol_));
    //if(!refine) cout << "No refinement level " << level << " at " << c.getcenter() << c.getsize() << endl;
    return refine;
  };
  virtual bool unify_cell_p (T_grid::const_local_iterator begin,
			     T_grid::const_local_iterator end) {
    return false;
  };
  virtual T_data get_data (const T_cell& c) {
    return T_data(get_rho(c), vec3d(0,0,0), T_lambda(n_lambda));
  }
    
  virtual int n_threads () {return 1;};
  virtual int work_chunk_levels () {return 5;};
  virtual int estimated_levels_needed () {return 0;};
};

void pascucci_test (po::variables_map opt)
{
  const int n_threads = opt["n_threads"].as<int>();
  const long n_rays =opt["n_rays"].as<long>();
  const int n_rays_temp = opt["n_rays"].as<long>();
  const vec3i grid_n (20,20,40);
  const int maxlevel = opt["maxlevel"].as<int>();
  const vec3d grid_extent(1000,1000,1000);
  const int n_forced=opt["n_forced"].as<int>();
  const int n_ir_forced=opt["n_ir_forced"].as<int>();
  const T_float I_rr=opt["i_rr"].as<T_float>();
  const int vis_ref_lambda = opt["vis_ref"].as<int>();
  const int ir_ref_lambda = opt["ir_ref"].as<int>();
  const T_float cameradist = 1e5;
  const string output_file(opt["output-file"].as<string>());
  const int seed=opt["seed"].as<int>();

  // ir equilibrium tolerances
  const T_float abs_tol = opt["abs_tol"].as<T_float>();
  const T_float rel_tol = opt["rel_tol"].as<T_float>();

  // opacity for the WD01 dust = 2.41e-6 kpc^2/Msun (=1.025e11 au^2/Msun)
  // for pascucci grid, n=rho0*rd*6.9078
  // which implies tau = 3.54e14 au^3/Msun * rho0

  // grid refinement tolerances
  const T_float r_tol =opt["r_tol"].as<T_float>();
  const T_float tau = opt["tau"].as<T_float>();
  const T_float tau_tol = opt["tau_tol"].as<T_float>();

  mcrx::seed(seed);
  T_unit_map units;
  units ["length"] = "au";
  units ["mass"] = "Msun";
  units ["luminosity"] = "W";
  units ["wavelength"] = "m";
  units ["L_lambda"] = "W/m";

  // read wavelengths
  //ifstream lfile("pascucci.wavelengths");
  ifstream lfile("pascucci_v2.wavelengths");
  vector<T_float> lambdas((istream_iterator<float>(lfile)),
			istream_iterator<float>());
  array_1 lambda(&lambdas[0], TinyVector<int,1>(lambdas.size()), neverDeleteData);
  const int lambda550 = lower_bound(lambdas.begin(), lambdas.end(),
				    0.551e-6) - lambdas.begin()-1;
  cout << "550nm is lambda entry " << lambda550 << endl;

  //need to communicate number of threads to the grain objects.
  Preferences prefs;
  prefs.setValue("n_threads", n_threads);
  prefs.setValue("use_grain_temp_lookup", true);

  // load grain model
  T_dust_model::T_scatterer_vector sv;
  sv.push_back(boost::shared_ptr<T_scatterer> 
	       //	       (new wd01_grain_model<polychromatic_scatterer_policy, mcrx_rng_policy>("/home/patrik/dust_data/crosssections", .1e-6, units)));
	       (new p04_grain_model<polychromatic_scatterer_policy, mcrx_rng_policy>("/home/patrik/dust_data/crosssections", prefs, units)));
  T_dust_model model (sv);
  model.set_wavelength(lambda); // this also resamples the grain_model objects
  // Convert optical depth to density
  const T_float unit_opacity_rho = 1.0/sv[0]->opacity()(lambda550);
  T_densities rho(1);
  rho = tau*unit_opacity_rho/(500*6.9078);
  T_densities n_tol(1); n_tol = tau_tol*unit_opacity_rho;
  cout << "Unit opacity density: " << unit_opacity_rho <<  ", n_tol: " << n_tol(0) << endl;

  // build grid
  cout << "Building grid" << endl;
  pascucci_factory factory (rho,n_tol, r_tol, maxlevel, 
			    opt["do_ir"].as<bool>() ? lambda.size() : 0);

  // first create the adaptive grid
  boost::shared_ptr<adaptive_grid<T_grid::T_data> > gg
    (new adaptive_grid<T_grid::T_data> (-grid_extent, grid_extent,
				grid_n, factory));

  // and then the full_sed_grid
  T_grid gr(gg, units);

  // calculate total dust mass
  T_densities mtot(1);
  mtot=0;
  for(T_grid::const_iterator c=gr.begin(); c!=gr.end(); ++c)
    mtot+=c->data()->densities()*c->volume();
  cout << "Total dust mass: " << mtot(0) << endl;

  cout << "\nSetting up emission" << endl;
  blackbody bb(5800,3.91e26);
  typedef emission<polychromatic_policy, T_rng_policy> T_emission;
  auto_ptr<T_emission> em;
  em.reset(new pointsource_emission<polychromatic_policy, T_rng_policy> (vec3d (.01, .01, -.01), bb.emission(lambda)));
  
  cout << bb.emission(lambda); 
  cout << "setting up emergence" << endl;
  std::vector<std::pair<T_float, T_float> > cam_pos;
  cam_pos.push_back(std:: make_pair (12.5*constants::pi/180, 0.1) );
  cam_pos.push_back(std:: make_pair (42.5*constants::pi/180, 0.1) );
  cam_pos.push_back(std:: make_pair (77.5*constants::pi/180, 0.1) );
  //cam_pos.push_back(std:: make_pair (78.5*constants::pi/180, 0.1) );
  //cam_pos.push_back(std:: make_pair (76.5*constants::pi/180, 0.1) );
  //cam_pos.push_back(std:: make_pair (77.5*constants::pi/180, .393) );
  //cam_pos.push_back(std:: make_pair (77.5*constants::pi/180, .785) );
  auto_ptr<T_emergence> e(new T_emergence(cam_pos, cameradist,
					  2.1*grid_extent[0], 300));

  cout << "shooting" << endl;

  std::vector<T_rng::T_state> states= generate_random_states (n_threads);
  T_rng_policy rng;

  //xfer<T_dust_model, polychromatic_grid, T_rng_policy> x (gr, *em, *e, model, rng, I_rr, 0,n_forced,10.,false, opt["do_ir"].as<bool>());
  xfer<T_dust_model, T_grid, T_rng_policy> x (gr, *em, *e, model, rng, I_rr, 0,n_forced,10.,false, opt["do_ir"].as<bool>());
  long n_rays_actual = 0;
  dummy_terminator t;
  T_biaser b(vis_ref_lambda);
  shoot (x, scatter_shooter<T_biaser> (b), t, states, n_rays,
  	 n_rays_actual, n_threads,0, "");
  T_float normalization = 1./n_rays_actual;
  
  {
    FITS output("!"+output_file, Write);
    e->write_images(output,normalization,units, "-STAR", false, false);

    // calculate integrated SEDs
    Table* iq_hdu=output.addTable("INTEGRATED_QUANTITIES",0);
    iq_hdu->addColumn(Tdouble, "lambda", 1,
		      units.get("L_lambda"));
    write(iq_hdu->column("lambda"),lambda,1);
    for (int i=0; i<cam_pos.size(); ++i) {
      integrate_image(output,
		      "CAMERA"+lexical_cast<string>(i)+"-STAR",
		      4*constants::pi*cameradist*cameradist,
		      "L_lambda"+lexical_cast<string>(i),
		      lambda,
		      "",
		      units,
		      dummy_terminator());
    }
  }

  cout << "Dust_model absorbed stellar luminosity: " 
       << sum(model.absorbed()*delta_quantity(lambda))*normalization << endl;

  if(!opt["do_ir"].as<bool>())
    return;

  // now we start the interesting stuff...
  // make dust mass vector
  array_2 dust_mass(gr.n_cells(), model.n_scatterers());
  array_2 star_intensity(gr.intensities().shape());

  int cc=0;
  const T_float area_factor = gr.area_factor();
  for(T_grid::const_iterator c=gr.begin(); c!=gr.end(); ++c, ++cc) {
    dust_mass(cc, Range::all()) = c->data()->densities()*c->volume();
    star_intensity(cc, Range::all()) =
      (gr.intensities()(cc,Range::all())*normalization/
       (4*constants::pi*c->volume()*area_factor));
  }

  //pascucci_factory ir_factory  (rho,n_tol, r_tol, maxlevel, lambda.size());

  // first create the adaptive grid
  boost::shared_ptr<adaptive_grid<T_ir_grid::T_data> > irgg
    (new adaptive_grid<T_ir_grid::T_data> (-grid_extent, grid_extent,
					   grid_n, gr.get_structure()));

  // and then the full_sed_grid
  T_ir_grid irg(irgg, dust_mass, lambda, units);

  array_2 dust_intensity(gr.n_cells(), lambda.size());
  dust_intensity=0;

  T_biaser b2(ir_ref_lambda);
  scatter_shooter<T_biaser> s(b2);
 
  // Determine dust equilibrium
  if (opt["do_equilibrium"].as<bool>())
    determine_dust_equilibrium_intensities
      (model,
       gr,
       irg,
       lambda,
       star_intensity,
       dust_intensity,
       states,
       n_rays_temp,
       n_threads,
       rel_tol,
       abs_tol,
       I_rr,
       n_ir_forced,
       10,
       s,
       t,
       true);
 
  // now generate output
  std::vector<grain_model<polychromatic_scatterer_policy, 
    mcrx_rng_policy>*> grain_models;
  for(int i=0; i< model.n_scatterers(); ++i) {
    grain_models.push_back
      (dynamic_cast<grain_model<polychromatic_scatterer_policy, 
       mcrx_rng_policy>*>(&model.get_scatterer(i)));
  }

  e.reset(new T_emergence(cam_pos, cameradist, 2.1*grid_extent[0], 50));
  irg.calculate_SED(array_2(star_intensity + dust_intensity), grain_models,
		    t, 1);//n_threads);
  gr.reset();
  xfer<T_dust_model, T_grid, T_rng_policy> x2 (gr, irg, *e, model, rng, 1e-2, 0,0,10.,false,false);

  cout << "Shooting dust emission" << endl;

  n_rays_actual = 0;
  shoot (x2, s, t, states, n_rays,
  	 n_rays_actual, n_threads, 0, "");

  normalization = irg.total_emission_weight()/n_rays_actual;
  FITS output(output_file, Write);
  e->write_images(output,normalization,units, "-IR", false, false);

  // calculate integrated SEDs
  for (int i=0; i<cam_pos.size(); ++i) {
    integrate_image(output,
		    "CAMERA"+lexical_cast<string>(i)+"-IR",
		    4*constants::pi*cameradist*cameradist,
		    "L_lambda_ir"+lexical_cast<string>(i),
		    lambda,
		    "",
		    units,
		    dummy_terminator());
  }
  
  /*
  {
    ofstream out("cell_seds.txt");
    for(T_ir_grid::iterator c=irg.begin(); c!=irg.end(); ++c)
      out << sqrt(dot(c->getcenter(),c->getcenter())) << '\n';
    out << "\n\n";
    for(T_ir_grid::iterator c=irg.begin(); c!=irg.end(); ++c) {
      for(int l=0; l<lambda.size(); ++l)
	out << lambda(l) << '\t' << c->data()->get_emission()(l) << '\n';
      out << "\n\n";
    }
  }
  */
  {
    ofstream out("gridstruct.txt");
    for(T_ir_grid::iterator c=irg.begin(); c!=irg.end(); ++c)
      out << c->getcenter()[0] << '\t' << c->getcenter()[1] << '\t'
	  << c->getcenter()[2] << '\n';
  }
}


int main (int argc, char** argv)
{
  // Declare the supported options.
  po::options_description desc("Pascucci et al 04 benchmark. Allowed options");
  desc.add_options()
    ("output-file", po::value<string>(), "output file name")
    ("n_rays", po::value<long>()->default_value(10000000), "number of rays")
    ("n_threads", po::value<int>()->default_value(8), "number of threads")
    ("maxlevel", po::value<int>()->default_value(13), 
     "max grid refinement level")
    ("n_forced", po::value<int>()->default_value(1), 
     "number of forced scatterings")
    ("n_ir_forced", po::value<int>()->default_value(0), 
     "number of forced scatterings for IR emission")
    ("vis_ref", po::value<int>()->default_value(17),
     "stellar radiation reference wavelength")
    ("ir_ref", po::value<int>()->default_value(46),
     "dust emission reference wavelength")
    ("i_rr", po::value<T_float>()->default_value(1e-2),
     "intensity where Russian Roulette starts")
    ("abs_tol", po::value<T_float>()->default_value(1e-6), 
     "absolute ir equilibrium tolerance")
    ("rel_tol", po::value<T_float>()->default_value(1e-2),
     "relative ir equilibrium tolerance")
    ("r_tol", po::value<T_float>()->default_value(0.2),
     "cell radial extent tolerance")
    ("tau", po::value<T_float>()->default_value(100),
     "edge-on V-band optical depth of disk")
    ("tau_tol", po::value<T_float>()->default_value(1),
     "cell optical depth refinement tolerance")
    ("seed", po::value<int>()->default_value(42), "random number seed")
    ("do_ir", po::value<bool>()->default_value(true), "calculate IR emission")
    ("do_equilibrium", po::value<bool>()->default_value(true),
     "calculate dust self-absorption")
    ("help", "print this page")
    ;

  po::positional_options_description p;
  p.add("output-file", 1);

  po::variables_map opt;
  po::store(po::command_line_parser(argc, argv).
          options(desc).positional(p).run(), opt);
  po::notify(opt);    

  if (opt.count("help")) {
    cout << desc << "\n";
    return 1;
  }
  if(!opt.count("output-file")) {
    cout << "Must specify output file name on command line" << endl;
    return 1;
  }

  cout << "Running Pascucci et al 04 benchmark with parameters:\n";
  cout << "output-file = " << opt["output-file"].as<string>() << '\n'
       << "n_rays = " << opt["n_rays"].as<long>() << '\n'
       << "n_threads = " << opt["n_threads"].as<int>() << '\n'
       << "maxlevel = " << opt["maxlevel"].as<int>() << '\n'
       << "n_forced = " << opt["n_forced"].as<int>() << '\n'
       << "n_ir_forced = " << opt["n_ir_forced"].as<int>() << '\n'
       << "vis_ref = " << opt["vis_ref"].as<int>() << '\n'
       << "ir_ref = " << opt["ir_ref"].as<int>() << '\n'
       << "i_rr = " << opt["i_rr"].as<T_float>() << '\n'
       << "abs_tol = " << opt["abs_tol"].as<T_float>() << '\n'
       << "rel_tol = " << opt["rel_tol"].as<T_float>() << '\n'
       << "r_tol = " << opt["r_tol"].as<T_float>() << '\n'
       << "i_rr = " << opt["i_rr"].as<T_float>() << '\n'
       << "tau = " << opt["tau"].as<T_float>() << '\n'
       << "tau_tol = " << opt["tau_tol"].as<T_float>() << '\n'
       << "seed = " << opt["seed"].as<int>() << '\n'
       << "do_ir = " << opt["do_ir"].as<bool>() << '\n'
       << "do_equilibrium = " << opt["do_equilibrium"].as<bool>() << "\n\n";

  pascucci_test ( opt );
}
